Code
import geopandas as gpd

#PASDA libraries shapefile
libraries = gpd.read_file("/Users/macbook/Desktop/Python/Week 14/PhiladelphiaLibraries201302/PhiladelphiaLibraries201302.shp")

#Swiming Pools 
pools = gpd.read_file("/Users/macbook/Desktop/Python/Week 14/PPR_Swimming_Pools/PPR_Swimming_Pools.shp")

#Playgrounds and Rec Centers
playgrounds = gpd.read_file("/Users/macbook/Desktop/Python/Week 14/PPR_Playgrounds/PPR_Playgrounds.shp")

#Parks
parks = gpd.read_file("/Users/macbook/Desktop/Python/Week 14/PPR_Program_Sites.geojson") 
Code
for name, gdf in {
    "libraries": libraries,
    "pools": pools,
    "playgrounds": playgrounds,
    "parks": parks
}.items():
    print(name, gdf.crs)
libraries EPSG:2272
pools EPSG:3857
playgrounds EPSG:3857
parks EPSG:4326
Code
TARGET_CRS = "EPSG:2272"

layers = {
    "libraries": libraries,
    "pools": pools,
    "playgrounds": playgrounds,
    "parks": parks
}

for name, gdf in layers.items():
    if gdf.crs != TARGET_CRS:
        layers[name] = gdf.to_crs(TARGET_CRS)

libraries = layers["libraries"]
pools = layers["pools"]
playgrounds = layers["playgrounds"]
parks = layers["parks"]
Code
libraries["type"] = "Library"
pools["type"] = "Pool"
playgrounds["type"] = "Playground"
parks["type"] = "Park"
Code
import pandas as pd
import geopandas as gpd
import folium
import censusdata
Code
u18_vars = [
    "B01001_001E",  # total population
    "B01001_003E","B01001_004E","B01001_005E","B01001_006E",
    "B01001_027E","B01001_028E","B01001_029E","B01001_030E",
]

acs_u18 = censusdata.download(
    "acs5",
    2022,
    censusdata.censusgeo([("state","42"), ("county","101"), ("tract","*")]),
    u18_vars
).reset_index()
Code
edu_vars = [
    "B15003_001E",  # total 25+
    "B15003_022E","B15003_023E","B15003_024E","B15003_025E"
]

acs_edu = censusdata.download(
    "acs5",
    2022,
    censusdata.censusgeo([("state","42"), ("county","101"), ("tract","*")]),
    edu_vars
).reset_index()
Code
inc_vars = ["B19013_001E"]

acs_inc = censusdata.download(
    "acs5",
    2022,
    censusdata.censusgeo([("state","42"), ("county","101"), ("tract","*")]),
    inc_vars
).reset_index()
Code
def make_geoid(df):
    df["GEOID"] = df["index"].apply(
        lambda x: "".join([part[1] for part in x.geo])
    )
    return df
Code
acs_u18 = make_geoid(acs_u18)
acs_edu = make_geoid(acs_edu)
acs_inc = make_geoid(acs_inc)
Code
acs_u18["under_18"] = acs_u18[
    [
        "B01001_003E","B01001_004E","B01001_005E","B01001_006E",
        "B01001_027E","B01001_028E","B01001_029E","B01001_030E"
    ]
].sum(axis=1)

acs_u18["pct_under_18"] = (
    acs_u18["under_18"] / acs_u18["B01001_001E"] * 100
)
Code
acs_edu["bachelors_plus"] = acs_edu[
    ["B15003_022E","B15003_023E","B15003_024E","B15003_025E"]
].sum(axis=1)

acs_edu["pct_bachelors_plus"] = (
    acs_edu["bachelors_plus"] / acs_edu["B15003_001E"] * 100
)
Code
acs_inc = acs_inc.rename(columns={
    "B19013_001E": "median_income"
})
Code
acs_full = (
    acs_u18[["GEOID", "pct_under_18"]]
    .merge(
        acs_edu[["GEOID", "pct_bachelors_plus"]],
        on="GEOID",
        how="left"
    )
    .merge(
        acs_inc[["GEOID", "median_income"]],
        on="GEOID",
        how="left"
    )
)
Code
import geopandas as gpd

# Read tracts geometry and merge the consolidated ACS table (acs_full)
tracts_geom = gpd.read_file(
    "https://www2.census.gov/geo/tiger/TIGER2022/TRACT/tl_2022_42_tract.zip"
)

tracts_geom = tracts_geom[tracts_geom["COUNTYFP"] == "101"]

tracts = tracts_geom.merge(
    acs_full,
    on="GEOID",
    how="left",
)
Code
# --- social infrastructure attribute renames and combine ---
libraries = libraries.rename(columns={"ASSET_NAME": "name"})
pools = pools.rename(columns={"amenity_na": "name"})
playgrounds = playgrounds.rename(columns={"park_name": "name"})
parks = parks.rename(columns={"park_name": "name"})
Code
social_infra = gpd.GeoDataFrame(
    pd.concat([
        libraries[["name", "type", "geometry"]],
        pools[["name", "type", "geometry"]],
        playgrounds[["name", "type", "geometry"]],
        parks[["name", "type", "geometry"]],
    ], ignore_index=True),
    crs=libraries.crs
)
Code
social_infra_web = social_infra.to_crs("EPSG:4326")
Code
# Project tracts and points to a projected CRS
tracts_proj = tracts.to_crs("EPSG:2272")
social_infra_proj = social_infra.to_crs("EPSG:2272")
Code
infra_join = gpd.sjoin(
    social_infra_proj,
    tracts_proj[["GEOID", "geometry"]],
    how="left",
    predicate="within"
)
Code
infra_counts = (
    infra_join
    .groupby("GEOID")
    .size()
    .reset_index(name="facility_count")
)
Code
infra_counts.head()
infra_counts["facility_count"].describe()
count    225.000000
mean       3.382222
std        2.032250
min        1.000000
25%        2.000000
50%        3.000000
75%        4.000000
max       13.000000
Name: facility_count, dtype: float64
Code
type_cols = {
    "Library": "library_count",
    "Playground": "rec_center_count",
    "Rec Center": "rec_center_count",
    "Park": "park_count",
    "Pool": "pool_count"
}
infra_by_type = (
    infra_join
    .groupby(["GEOID", "type"])
    .size()
    .unstack(fill_value=0)
    .reset_index()
    .rename(columns=type_cols)
)
for col in type_cols.values():
    if col not in infra_by_type.columns:
        infra_by_type[col] = 0
Code
tracts_access = (
    tracts_proj
    .merge(infra_counts, on="GEOID", how="left")
    .merge(infra_by_type, on="GEOID", how="left")
)

# Replace NA with 0
count_cols = ["facility_count", "library_count", "rec_center_count", "park_count", "pool_count"]
tracts_access[count_cols] = tracts_access[count_cols].fillna(0)
Code
tracts_access_web = tracts_access.to_crs("EPSG:4326")
Code
count_field_options = [
    ("Social Infrastructure Count:", ["facility_count"]),
    ("Libraries:", ["library_count", "Library"]),
    ("Rec Centers:", ["rec_center_count", "Rec Center", "Playground"]),
    ("Parks:", ["park_count", "Park"]),
    ("Pools:", ["pool_count", "Pool"]),
]
tooltip_fields = []
tooltip_aliases = []
for alias, options in count_field_options:
    for col in options:
        if col in tracts_access_web.columns:
            tooltip_fields.append(col)
            tooltip_aliases.append(alias)
            break
if not tooltip_fields:
    tooltip_fields = ["facility_count"]
    tooltip_aliases = ["Social Infrastructure Count:"]
Code
from folium.plugins import DualMap

dual_map = DualMap(
    location=[39.9526, -75.1652],
    zoom_start=11,
    tiles="CartoDB positron"
)
Code
import numpy as np

def style_facilities(feature):
    v = feature["properties"].get("facility_count")
    if v is None or v == 0 or (isinstance(v, float) and np.isnan(v)):
        return {
            "fillColor": "#f0f0f0",
            "color": "#cccccc",
            "weight": 0.3,
            "fillOpacity": 0.45
        }
    return {
        "fillColor": cmap_facilities(v),
        "color": "#ffffff",
        "weight": 0.3,
        "fillOpacity": 0.9
    }
Code
tracts_access_web["facility_count"].describe()
count    408.000000
mean       1.865196
std        2.260431
min        0.000000
25%        0.000000
50%        1.000000
75%        3.000000
max       13.000000
Name: facility_count, dtype: float64
Code
import branca.colormap as cm


bins = tracts_access_web["facility_count"].quantile(
    [0, 0.2, 0.4, 0.6, 0.8, 1]
).values

cmap_facilities = cm.StepColormap(
    colors=cm.linear.Blues_06.colors,
    index=bins,
    vmin=bins[0],
    vmax=bins[-1]
)
Code
import censusdata
import geopandas as gpd
import pandas as pd
import folium
Code
acs_vars = {
    "total_pop": "B01001_001E",
    "male_under_18": [
        "B01001_003E","B01001_004E","B01001_005E","B01001_006E"
    ],
    "female_under_18": [
        "B01001_027E","B01001_028E","B01001_029E","B01001_030E"
    ]
}

tracts_acs = censusdata.download(
    "acs5",
    2022,
    censusdata.censusgeo([("state","42"), ("county","101"), ("tract","*")]),
    [acs_vars["total_pop"]] +
    acs_vars["male_under_18"] +
    acs_vars["female_under_18"]
)
Code
tracts_acs["under_18"] = tracts_acs[
    [
        "B01001_003E",
        "B01001_004E",
        "B01001_005E",
        "B01001_006E",
        "B01001_027E",
        "B01001_028E",
        "B01001_029E",
        "B01001_030E",
    ]
].sum(axis=1)
Code
edu_vars = [
    "B15003_001E",  # total 25+
    "B15003_022E",  # bachelor's
    "B15003_023E",  # master's
    "B15003_024E",  # professional
    "B15003_025E"   # doctorate
]
Code
income_vars = ["B19013_001E"]
Code
tracts_acs_edu = censusdata.download(
    "acs5",
    2022,
    censusdata.censusgeo([("state","42"), ("county","101"), ("tract","*")]),
    edu_vars
)

tracts_acs_inc = censusdata.download(
    "acs5",
    2022,
    censusdata.censusgeo([("state","42"), ("county","101"), ("tract","*")]),
    income_vars
)
Code
tracts_acs_edu = tracts_acs_edu.reset_index()
tracts_acs_inc = tracts_acs_inc.reset_index()
Code
tracts_acs_edu = tracts_acs_edu.rename(columns={
    "B15003_001E": "total_25plus",
    "B15003_022E": "bachelors",
    "B15003_023E": "masters",
    "B15003_024E": "professional",
    "B15003_025E": "doctorate"
})
Code
# If total population is still B01001_001E, rename it first
tracts_acs = tracts_acs.rename(columns={
    "B01001_001E": "total_pop"
})

# Recreate under_18 if needed
if "under_18" not in tracts_acs.columns:
    tracts_acs["under_18"] = tracts_acs[
        [
            "B01001_003E","B01001_004E","B01001_005E","B01001_006E",
            "B01001_027E","B01001_028E","B01001_029E","B01001_030E"
        ]
    ].sum(axis=1)

# NOW create the percentage
tracts_acs["pct_under_18"] = (
    tracts_acs["under_18"] / tracts_acs["total_pop"] * 100
)
Code

if "bachelors" in tracts_acs_edu.columns and "pct_bachelors" not in tracts_acs_edu.columns:
    tracts_acs_edu["pct_bachelors"] = (
        tracts_acs_edu["bachelors"]
        / tracts_acs_edu["total_25plus"] * 100
    )
Code
tracts_acs_inc = tracts_acs_inc.rename(columns={
    "B19013_001E": "median_income"
})
Code
tracts_acs_edu["GEOID"] = tracts_acs_edu["index"].apply(
    lambda x: "".join([part[1] for part in x.geo])
)
Code
tracts_acs_inc["GEOID"] = tracts_acs_inc["index"].apply(
    lambda x: "".join([part[1] for part in x.geo])
)
Code
tracts_acs["GEOID"] = tracts_acs.index.map(
    lambda x: (
        x.geo[0][1] +  # state
        x.geo[1][1] +  # county
        x.geo[2][1]    # tract
    )
)
Code
tracts_acs[["total_pop", "under_18", "pct_under_18"]].head()
total_pop under_18 pct_under_18
Census Tract 1.01; Philadelphia County; Pennsylvania: Summary level: 140, state:42> county:101> tract:000101 1947 40 2.054443
Census Tract 1.02; Philadelphia County; Pennsylvania: Summary level: 140, state:42> county:101> tract:000102 2897 162 5.591992
Census Tract 2; Philadelphia County; Pennsylvania: Summary level: 140, state:42> county:101> tract:000200 3486 414 11.876076
Census Tract 3; Philadelphia County; Pennsylvania: Summary level: 140, state:42> county:101> tract:000300 3914 280 7.153807
Census Tract 4.01; Philadelphia County; Pennsylvania: Summary level: 140, state:42> county:101> tract:000401 2675 127 4.747664
Code
acs_full = (
    tracts_acs[["GEOID", "pct_under_18"]]
    .merge(
        tracts_acs_edu[["GEOID", "pct_bachelors"]],
        on="GEOID",
        how="left"
    )
    .merge(
        tracts_acs_inc[["GEOID", "median_income"]],
        on="GEOID",
        how="left"
    )
)
Code
tracts = tracts_geom.merge(
    acs_full,
    on="GEOID",
    how="left"
)

tracts_web = tracts.to_crs("EPSG:4326")
Code
import folium

m_acs = folium.Map(
    location=[39.9526, -75.1652],
    zoom_start=11,
    tiles="CartoDB positron"
)
Code
import numpy as np

bins_income = tracts_web["median_income"].quantile(
    [0, 0.2, 0.4, 0.6, 0.8, 1]
).values
Code
import branca.colormap as cm

cmap_income = cm.StepColormap(
    colors=cm.linear.YlOrRd_06.colors,
    vmin=bins_income[0],
    vmax=bins_income[-1],
    index=bins_income
)
Code
bins_bach = tracts_access_web["pct_bachelors_plus"].quantile(
    [0, .2, .4, .6, .8, 1]
).values
Code
import branca.colormap as cm
import numpy as np

cmap_bach = cm.StepColormap(
    colors=cm.linear.Purples_06.colors,
    index=bins_bach,
    vmin=bins_bach[0],
    vmax=bins_bach[-1]
)
cmap_bach.caption = "% Bachelor's Degree or Higher"
Code
def style_bachelors(feature):
    v = feature["properties"].get("pct_bachelors_plus")
    if v is None or np.isnan(v):
        return {
            "fillColor": "#A9A9A9",
            "color": "white",
            "weight": 0.4,
            "fillOpacity": 0.85
        }
    return {
        "fillColor": cmap_bach(v),
        "color": "white",
        "weight": 0.4,
        "fillOpacity": 0.85
    }
Code
folium.GeoJson(
    tracts_web,
    name="% Bachelor’s Degree or Higher",
    style_function=style_bachelors
).add_to(m_acs)
<folium.features.GeoJson at 0x14ee4e5d0>
Code
cmap_income = cm.linear.YlOrRd_09.scale(
    tracts_web["median_income"].min(),
    tracts_web["median_income"].max()
)
Code
import numpy as np

bins_income = tracts_web["median_income"].quantile(
    [0, 0.2, 0.4, 0.6, 0.8, 1]
).values
Code
import branca.colormap as cm

cmap_income = cm.StepColormap(
    colors=cm.linear.YlOrRd_06.colors,
    vmin=bins_income[0],
    vmax=bins_income[-1],
    index=bins_income
)
Code
def build_acs_layer(name, field, palette):
    values = tracts_access_web[field]
    bins = values.quantile([0, 0.2, 0.4, 0.6, 0.8, 1]).values
    cmap = cm.StepColormap(
        colors=palette,
        index=bins,
        vmin=bins[0],
        vmax=bins[-1]
    )
    cmap.caption = name
    def style_fn(feature, field=field, cmap=cmap):
        val = feature["properties"].get(field)
        if val is None or (isinstance(val, float) and np.isnan(val)):
            return {"fillColor": "#d9d9d9", "color": "white", "weight": 0.3, "fillOpacity": 0.4}
        return {"fillColor": cmap(val), "color": "white", "weight": 0.3, "fillOpacity": 0.85}
    return {"name": name, "field": field, "style": style_fn, "cmap": cmap}

acs_layers = [
    build_acs_layer("% Bachelor's Degree or Higher", "pct_bachelors_plus", cm.linear.Purples_06.colors),
    build_acs_layer("% Under Age 18", "pct_under_18", cm.linear.Greens_06.colors),
    build_acs_layer("Median Household Income", "median_income", cm.linear.YlOrRd_06.colors),
]

FINAL MAP

Code
import folium
from folium.plugins import DualMap
from folium.features import GeoJsonTooltip
import branca.colormap as cm
import numpy as np
Code
#Dual Map Setup

from folium.plugins import DualMap
import folium
from folium.features import GeoJsonTooltip

dual_map_edu = DualMap(
    location=[39.9526, -75.1652],
    zoom_start=11,
    tiles="CartoDB positron"
)
Code
folium.GeoJson(
    tracts_access_web,
    style_function=style_facilities,
    tooltip=GeoJsonTooltip(
        fields=tooltip_fields,
        aliases=tooltip_aliases,
        localize=True,
        sticky=True
    )
).add_to(dual_map_edu.m1)

# Legend (optional but fine here)
cmap_facilities.add_to(dual_map_edu.m1)
0.02.24.36.58.710.813.0
Code
for layer in acs_layers:
    fg = folium.FeatureGroup(name=layer["name"], show=(layer["field"] == "pct_bachelors_plus"))
    folium.GeoJson(
        tracts_access_web,
        style_function=layer["style"],
        tooltip=GeoJsonTooltip(
            fields=[layer["field"]],
            aliases=[f"{layer['name']}:"],
            localize=True,
            sticky=True
        )
    ).add_to(fg)
    fg.add_to(dual_map_edu.m2)
    layer["cmap"].add_to(dual_map_edu.m2)

folium.LayerControl(collapsed=False).add_to(dual_map_edu.m2)
0.529941706412294716.532.548.464.480.396.29629629629629% Bachelor's Degree or Higher
Code
tracts_access_web.head()
STATEFP COUNTYFP TRACTCE GEOID NAME NAMELSAD MTFCC FUNCSTAT ALAND AWATER ... INTPTLON geometry pct_under_18 pct_bachelors_plus median_income facility_count Library Park Playground Pool
0 42 101 026301 42101026301 263.01 Census Tract 263.01 G5020 S 406340 0 ... -075.1637161 POLYGON ((-75.16899 40.07147, -75.16864 40.071... 22.435737 28.359827 57992 0.0 0.0 0.0 0.0 0.0
1 42 101 029200 42101029200 292 Census Tract 292 G5020 S 1745920 28444 ... -075.1017831 POLYGON ((-75.11288 40.02746, -75.1128 40.0274... 23.318995 19.735431 82784 0.0 0.0 0.0 0.0 0.0
2 42 101 024400 42101024400 244 Census Tract 244 G5020 S 421189 0 ... -075.1638925 POLYGON ((-75.16912 40.02387, -75.16843 40.024... 33.115265 40.127077 56912 3.0 0.0 1.0 2.0 0.0
3 42 101 033200 42101033200 332 Census Tract 332 G5020 S 842716 0 ... -075.0449684 POLYGON ((-75.05464 40.04465, -75.05396 40.045... 29.119272 26.290400 80409 2.0 0.0 1.0 0.0 1.0
4 42 101 980200 42101980200 9802 Census Tract 9802 G5020 S 4916641 250846 ... -075.0443487 POLYGON ((-75.07448 40.08685, -75.07443 40.087... 6.666667 24.867725 -666666666 5.0 0.0 1.0 4.0 0.0

5 rows × 21 columns

Code
dual_map_edu
Make this Notebook Trusted to load map: File -> Trust Notebook
Code
folium.GeoJson(
    tracts_access_web,
    style_function=style_facilities,
    tooltip=GeoJsonTooltip(
        fields=tooltip_fields,
        aliases=tooltip_aliases,
        localize=True,
        sticky=True
    )
).add_to(dual_map.m1)
<folium.features.GeoJson at 0x16ba98150>
Code
for layer in acs_layers:
    fg = folium.FeatureGroup(name=layer["name"], show=(layer["field"] == "pct_bachelors_plus"))
    folium.GeoJson(
        tracts_access_web,
        style_function=layer["style"],
        tooltip=GeoJsonTooltip(
            fields=[layer["field"]],
            aliases=[f"{layer['name']}:"],
            localize=True,
            sticky=True
        )
    ).add_to(fg)
    fg.add_to(dual_map.m2)
    layer["cmap"].add_to(dual_map.m2)

folium.LayerControl(collapsed=False).add_to(dual_map.m2)
<folium.features.GeoJson at 0x16b7a7d10>
Code
# Add aligned legends for both map panes
cmap_facilities.add_to(dual_map.m1)
0.529941706412294716.532.548.464.480.396.29629629629629% Bachelor's Degree or Higher
Code
dual_map
Make this Notebook Trusted to load map: File -> Trust Notebook